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ABSTRACT 

Long lived modes of elliptical galaxies can exist d la van Kampen. Specific systems 
may possess long lived oscillations which Landau damp on time scales longer than a 
Hubble time. Some physical processes such as a close encounter, tidal forces from a 
cluster or an orbiting satellite could preferentially excite a coherent mode. These may 
relate to the observed faint structure in elliptical galaxies such as shells and ripples. 
Their detection in projected phase space would ultimately provide a detailed probe of 
the underlying potential. 

I give an overview of linear perturbations to stationary solutions of the Vlasov 
equation, including a discretized Hermite polynomial expansion which explicitly 
demonstrates completeness and orthogonality of solutions. Some exact solutions are 
shown, which implies the feasibility of such a procedure and suggest future fully 
numerical studies. 



Subject headings: galaxies: internal motions, galaxies: structure, stars: stellar 
dynamics 



1. Introduction 



Elliptical galaxies provide a clean laboratory to study the physics of self-gravitating fluids. 
We would like to study the possibility of applying normal mode analysis to complement the 
observational data. The conventional wisdom is that most excitations Landau damp or phase 
mix at rates comparable to the dynamical time, making it impossible to observe the equivalent of 
p-modes in stars. Exceptions to this rule have been studied by previous authors. Mathur (1990) 
has shown that not all modes decay by proving the existence of discrete eigenvalues that allow 
stable oscillations. Weinberg (1991) numerically studied these modes in the vertical structure of 
the disk. A different approach was shown in (Weinberg 1993) where he numerically studied modes 
which decay on time scales long compared to the dynamical time. In this article I wish to extend 
this search by examining qualitative relationships between various modes in linear theory. 
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I follow the van Kampen approach to construct these modes. Since Case (1959) has shown 
van Kampen's method to be equivalent to Landau's treatment, we can use the qualitative features 
of van Kampen modes to infer properties of Landau damping. In the Jeans analysis, for example, 
the normal van Kampen modes are singular (Binney and Tremaine 1987, appendix 5A). This 
not only makes them impossible to excite, but also violates the Landau damping assumptions, 
where the distribution must be differentiable for v. Conversely, I will suggest how non-singular 
modes might be constructed, and show explicit examples in some one dimensional toy systems. A 
differentiable van Kampen mode implies the existence of a mode which does not Landau damp. If 
the system is furthermore proven to be stable, we know that the Landau integration contour must 
contain a pole on the real axis, without actually having to find it. 

In this paper I first give a general treatment of linear perturbations to review the notation and 
physical interpretation. Next I will assume spherical symmetry to study the modes found in simple 
stationary systems. While they are too simple to be applicable to real astrophysical systems, they 
demonstrate how non-singular van Kampen modes might exist. Following that, I will construct 
a frame work based on matrix expansions, which could be used to study modes numerically. 
This method is applied to the example of the Jeans' instability. We then construct an explicitly 
self-adjoint complete representation of the linear Vlasov operator in the presence of discontinuous 
functions. By showing how an unbounded operator which is not explicitly self-adjoint can possess 
a complete set of orthogonal eigenfunctions we find a resolution to Case's puzzle. The argument 
also lends strength to the feasibility of a general mode decomposition. The paper concludes with 
some speculation on possible astrophysical implications and suggestions for future work. 



2. Formulation 



Qualitatively, one proceeds as follows: we consider a stable stationary system (such as an 
elliptical galaxy), and look for periodic perturbations. This is accomplished by first constructing 
a small perturbation consisting of a smooth periodic pattern in the background potential, i.e. a 
periodic hamiltonian flow. To linear order, this pattern does not interact with itself. It does, 
however, effect the background fluid. Having assumed the system stable, we know that no resonant 
growth can result. We need only study the response of the background, and we are done. If 
we are unfortunate, the background response may cancel the perturbation exactly, and such an 
example is given in the next section. Or the background response to a periodic perturbation 
may be chaotic, for example if the potential is not integrable. In general, however, we have some 
non-trivial background response, as constructed explicitly in section ^ for the case of the Jean's 
instability. 
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To phrase this mathematically, we are interested in the solution to Vlasov's equation 



dt ^ 



1=1 



' dxi 



df 
dxi dvi 



(1) 
(2) 

The problem under consideration consists of a stationary background solution fh{x,v) and a 
periodic perturbation fp{x, v) such that f = fb + ee*'^*/p, where e is a small number. The sufficient 
condition for the perturbation equation to be valid is that 



V2$ = 47rGp, p{x)= fd^v, x = {xi,X2,X3), f = f{x,v;t). 



' dxi 



dxi 



(3) 



at all points x for i = 1,2, 3, which is weaker than requiring efp < fh at all points in phase space. 
The latter condition only becomes necessary when we consider negative perturbations, but positive 
perturbations are allowed even when singular in velocity space. 

To first order in e the perturbation satisfies the linear equation 



Li 



3 r 

i=l 



: Li + L2 



d 



' dxi dxi dvj 



i=l 



dvi dxj 



d:\. 



(4) 



This decomposition allows a simple physical interpretation. Solutions to Li describe the 
motion of the perturbation, and L2 describes the response of the background system. Mathur 
(1990) showed that for any oj which is not an eigenvalue of Li, L can inherit a discrete eigenvalue 
from L2. In this paper, however, we will consider the eigenvalues of Li. Any eigenvalue of Li 
corresponds to an eigenvalue of L since the velocity space structure of L2/p depends solely on the 
background distribution. To see this, let Li/^ = ojfa- We can construct the response fp such that 
the total perturbation fp = fa + fp- It satisfies the inhomogeneous equation 



(L - uj)fp 



dfb 



3 

y 

^ dvi dxi 



(5) 



Equation may be singular. We thus can not simply invert h — u>. Instead one must establish 
that at least one solution exists. This can be seen by going back to the linearized Vlasov 
equation (^), where we use the full time dependence to solve a linear hyperbolic problem in a 
periodically driven field 6*'^*$^: 



1=1 



dvi dxi 



(6) 
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Start with the initial value /^{t = 0) = 0, and we have a unique solution to the initial value 
problem. Wait for a time T, and Fourier transform the result, yielding f^{T,ujt). Repeat this for 
larger values of T, taking the limit as T — > oo. Since we are working with a stable system, there 
can be no growing modes with negative imaginary eigenvalues. Any damped modes with negative 
eigenvalues we simply discard. We take only the resonant frequency piece, and transform back to 
get fs = f-yiT = Qo,uJt = Lo)e^^^ . This will satisfy ^. To see this, take any interval [to, ^1=^0 + ^] 
with to > 0, and apply a discrete Fourier transform on that interval to equation ^. The RHS 
contains only one frequency component oj. Therefore the LHS operator = dt + Ij applied to the 
solution /e = Lfif^ will contain only that same frequency component. Since does not depend 
on t, it can only map functions to zero, but not generate new frequencies. Thus fs solves ^, and 

To numerically approximate (^), one can discretize the operators by expanding all functions 
in a given basis (Pen and Jiang 1992, Pen 1992). If we expand as a discrete sum along the lines 
of section we obtain a sequence of approximations which will converge to For each order 
in the expansion, the differential and integral operators are finite and discrete, so we can solve 
(||). We use the Gram-Schmidt orthogonalization to invert L — and require orthogonality to 
fa when we encounter the matrix singularity. Specific examples of exact and series solutions will 
follow below. 

We now only need to consider Li. Its solutions describe the motion of an ensemble of particles 
in a static field. A single periodic orbit corresponds to a localized distribution ((5-function) in phase 
space, which is a periodic solution, and thus a discrete sum of eigenstates. In the case of integrable 
systems, one can describe all eigenstates in terms of the action-angle coordinates. In general it is 
possible to have periodic patterns, whether or not individual orbits are periodic. The pattern need 
not have the same period as its constituent orbits. To illustrate this point, consider a single orbit 
in an integrable system with period T. A single particle corresponds to a distribution function 
with pattern period T, but n particles equispaced along the angle variable with have period T/n. 
If we fill the orbit with a continuum of particles, the pattern will be stationary. Now consider 
a second orbit with a different period. If the ratio of the two periods is a rational number, we 
can construct a pattern with period equal to any linear sum of the two constituent orbit periods. 
Given a set with a continuum parameter of orbital periods, we can, by judicious choice of orbits 
and phase angles, construct a pattern with any period we wish. In astrophysical situations, such 
patterns might be observable. We will now invert the procedure and try to construct patterns 
directly from the Eulerian description of the distribution function. 

To linear order in phase space, the Liouville theorem assures that these solutions do not 
diffuse or damp since dynamical friction is a higher order non-linear phenomenon. From (|3|) it 
follows that a constant mass density perturbation will have more linear behavior if it is smeared 
out in space. One would thus expect that long wavelength perturbations are dynamically longer 
lived, which will make their detection easier observations with limited spatial resolution. 
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3. Spherical Potentials 



In general, the eigenmodes of a system must be computed numericaUy. For a few systems, 
however, it is possible to find exact or series expressions of some modes. I will show some 
examples to illustrate the effects in simple models. They are radial toy models with non-singular 
van-Kampen modes, which implies the existence of initial conditions which do not Landau damp. 
While they do not have a clear correspondence to three dimensional real galaxies, they are easier 
to calculate and may provide some hints about the behaviour of realistic systems. 

Consider simple power law models, where the background distribution has a power law 
dependence on the radial coordinate p oc r". The simplest is n = which is realized in the 
Einstein sphere, see for example (Mikahilovskii 1972). The Einstein sphere is a stable distribution 
describing a constant density sphere where all particles are on tangential (non-radial) orbits, and 
the distribution function is isotropic in the tangent plane. All perturbation orbits which do not 
leave the sphere are periodic, and thus all modes have the same eigenvalue. This exemplifies the 
existence of an isolated eigenvalue in Li. 

For n = —1, the potential gradient V<l>b = F \s constant. Now consider purely radial 
perturbation orbits. The linear Vlasov equation ^ becomes 

-iUU = Vr^-F^^=iu;f,. (7) 

Laplace transforming fp with respect to r, such that fp = J fke~^^dk allows us to express the 
radial eigenmode as 

fk = eM=^+i^-k\r\)- (8) 

This solution is certainly smooth everywhere except at r = 0. We thus need to modify the 
background potential to allow a smooth transition. The 1/r density has rapidly divergent mass, 
so one must limit the power law at some radius R, which provides a characteristic scale for k. 
By using positive and negative perturbations of the lower harmonics, one can easily construct 
non-trivial waves and modes for standing radial density waves. 

The zero eigenvalues of isothermal spheres (or any other power law system) can also be solved. 
This scenario is attractive since many astrophysical objects have such a dependence. Let be the 
velocity dispersion of the background material. Consider an isothermal Maxwellian perturbation 
with some other velocity dispersion ap, as might describe an elliptical galaxy embedded in some 
halo. From Li/a = we obtain the well known result pa = r", n = —2ab/ap. The background 
(halo) response ph to the imbedding is 

2 

-I- on -I- 8 

and the net density fluctuation is Pp = Pa + Pb- In the case that ct), = ap, one would obtain n = —2. 
The background response is equal and opposite to the perturbation and exactly cancels it. In this 
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case, no net perturbation actually happened, except for relabeling b particles to belong to p. But 
if we solve the perturbation equation exactly, drV^h/p = becomes 

which implies n = (—5 it \/—7)/2, and we have a nontrivial solution. This eigenmode can be 
obtained from (^) in the limit as pr — > oo. 



4. Jeans Instability 



The Jeans instability is well understood with known analytic van Kampen modes. It is thus a 
good testbed for general analysis. The spatial translation invariance allows an exact series solution 
in terms of an infinite sequence of discrete matrix operators on successively refined basis functions. 

Consider a homogeneous isotropic Maxwellian background fluid whose distribution function 
is given as = ph ex.p{—v'^ /2a'^) / a\/^ . Then apply the Jeans swindle and set the potential of 
this component to zero. Let the Fourier mode fp be periodic in space with wave number k in one 
dimension (say x) and constant along the other two (y, z) so the perturbation equation reads 

U + ^kvU-'^^^e-^y^'|dvU = (11) 



where s = \f2a. Now L = fcf — (A;j/\/27r(7/i;)t; exp(— v^/s^) / dv with /cj = ArKGphjcy^ being the 
Jeans wavenumber. Equation ([Tl|) can be solved as 



U(t) = e^^Vp(O). (12) 



All we need to do is flnd the eigenvalues of L and project the initial condition /(O) onto the 
eigenstates to obtain the complete solution. 

The simplest case is the free equation, where G = so we have no gravity. Then the 
eigenmodes of L are /^q = d{v — vq), Dirac (^-functions which move with frequency uj = kvQ. The 
next simplest case are stationary solutions, i.e. distributions corresponding to zero eigenvalues, 
L/p = 0. Let fp = exp{ikx){a6{v) + bm{v)) where m{v) = exp(— 'L'^/s^)/s-^/7r is a Maxwellian 
distribution similar to fb- Equation ( |ll| ) becomes 

k \ a , , 

k]j -' = 1 <^^> 

and we have solved the static perturbation equation as being the sum of the Gaussian and a 6 
function. The relevant limits are: 
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• k ^ oo: for very short wavelength perturbations, 6^0 and our solutions are (5- functions, 
which solve the free equations as expected. 

• k ^ kj: at the Jeans' length, the solution is a plain Gaussian without any delta function, 
which is a static perturbation. 

• A: — > ^/2kj: this is the equality point, where half the density is in a J- function and half is in 
the Gaussian, so a = 6. 

• A; — > 0: it certainly seems curious that we have static solutions with wavelengths much longer 
than the Jeans' length. Note that a —b, which means that /9p — > and the net density 
fluctuation tends to zero. All the structure is in velocity space, implying that we can always 
have coherent stable perturbations, even though they contain less and less mass. 



We can now consider the complete spectral decomposition of L. From the static and free 
solutions we are led to the ansatz fp = ex.p{ikx){iJ,6{v — vq) + fh) and expand in Hermite 
polynomials 



i=0 



(14) 



where z/ = v/s, and we get scq = ph- The Hi = — 1,... are defined using the 

conventions of Gradshteyn and Ryzhik (1980), and the normalization Ni = {2^n\^/^T)^^^'^ . The 
completeness of this expansion for Lebesque integrable functions is shown in (Keener). Note 
that projection ( [l^ is onto a skew basis, where the inverse projection occurs through plain 
Hermite polynomials Ci = Ni J duHifh without an exponential weighting. 



First, set /x = 0. Substituting (14) into (|llD , multiplying by NiHi and integrating over du, we 
have turned the continuum equation (|ll]) into a discrete matrix equation, so writing the q as a 
column vector c, we obtain c = zLc where L now becomes 
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The eigenvalues are easily read off from L. If A; > kj, we can symmetrize entries L12 and L21 to 
their geometric mean (Wilkinson 1965), yielding all real eigenvalues, implying that all eigenmodes 
are stable. The apparent dispersion becomes clear: a density perturbation 6co projects onto all 
eigenstates, which propagate at different velocities. Only the norm of c is conserved, and cq clearly 
decays, explaining the dispersion in density space. For k = kj, the matrix becomes explicitly 
singular, with the Gaussian cq being a static eigenmode, as ana lyzed abov e. For k < kj , we obtain 
two (and only two) imaginary eigenvalues, uJi_2 ^ ±2iTT^^'^ka^k'^ /kj — 1. This result holds for 
small k and is obtained as leading order term from the continued fraction expansion described 
below, by directly observing the singularity of the matrix. There are only two unstable eigenstates, 
one growing and one decaying. 

Let us now construct eigenstates. The eigenstates of V are 5-functions, as is easily verified 
from the recurrence relation of Hermite polynomials. Unfortunately, convergence to such 
discontinuous functions is quite slow. We thus return to our ansatz and solve for the analytic part 
separately. Restoring fj, into our calculation, we can solve for the eigenstate iteratively. Substitute 
fi5{v — vq) back into equation (1TT|). We then obtain 



kvp \ ^_ vpkj 
r~ k 



( \ 






(16) 



which we can solve for c, unless A; = 0, which we had discussed earlier. Gaussian elimination is 
most easily applied from the bottom up, and we obtain for the top two coefficients 



Acq + ci =0 
S]e„+(A-C(A->))., = ^-j^ 



(17) 
(18) 



Here C(x) is defined from the continued fraction 
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while A = — fo/2(T7r^/^. Knowing the first two elements, the remainder are simply given by the 



tridiagonal recurrence relation (16). This algorithm can be verified by starting with a finite matrix 



and iteratively refining the solution. 

The same method can be applied to solve for the imaginary eigenvalue, where we simply 
require ( p!8| ) to be singular, and interpret A as an eigenvalue. We thus have an algebraically implicit 
expression for the growth factor. Since the imaginary eigenvalue is unique up to sign, a growing 
ansatz should obtain the correct solution and the implicit solution described here equivalent to 
equation (5-31) in Binney and Tremaine (1987). 

The qualitative features are easy to extract. As A: ^ oo or t>o ^ oo, the analytic component 
vanishes, c — > and we recover the free non-interacting solution. As /c ^ we obtain 
— Aci = Co = — /u unless we move vq to oo at the same time, in which case arbitrary ratios of cq/ci 
can be achieved. 



We can compare this to the analytic expression for the van Kampen mode 

6{v-uj/k) (20) 



~ (27ra2)3/2fc2 + 



where 

Wy(Z) = -^p / — /^cix (21) 



^/2tt J-oo X - Z 

and p denotes the Cauchy principal value for the integral. We verify correct convergence for 
k > kj. The matrix result, however, is general, and correctly yields all the continuum and discrete 
eigenvalues and eigenvectors for all values of k. 

This example illustrates the separation between the solutions of Li which are simply the 
6 functions, and the contribution from L2 which contains the response. This approach can be 
numerically extended by expanding the distribution function at every point in space and solving 
the more complex set of equations for arbitrary potentials. 



5. Discussion 



Formally, L appears self-adjoint, which unfortunately only holds for (differentiable) 
functions. But as we saw, the self-adjointness can also be elucidated for quite pathological function 
domains. In non-trivial background potentials, smoothness of eigenmodes should increase since Li 
contains a d^, which would diverge strongly for singular functions, such as occurred in the Jeans 
analysis. 

For general systems such as elliptical galaxies, there is no way to separate the distribution 
function into the product of a radial and a velocity piece, as we did in all the examples. This foils 
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any attempt to apply a Fourier or other integral transform to r and v separately to express the 
differential operators as algebraic ones. Already in the Jeans example, we saw that one cannot 
define dispersion measures for van Kampen modes. In the absence of new tools, the general 
problem is very difficult, which is to be expected of partial differential equations. Linearity allows 
systematic numerical sudies as we will see below. 

The assumption of being a small perturbation holds as long as the inequality is satisfied, 
which is even possible in the presence of a divergent (5-function. One need only make sure that the 
coefficient of that 5-function is positive to prevent any possibility of causing negative densities in 
phase space. 

Most elliptical galaxies are observed to possess some weak small scale structure, which should 
certainly be subject to perturbative modeling. If we can have a complete analytic or numerical 
understanding of the eigenmodes, these perturbations can teach us much about the detailed 
potential structure. A direct prediction of this analysis is that perturbations should posses a 
discrete symmetry: positive and negative density perturbations should occur with comparable 
frequencies, and display similar patterns. Linearity of (^) allows us to reverse the sign of the 
perturbation fp — > —fp, so we do not expect small perturbations to be skewed when averaged 
over many instantiations, i.e. < pp >= 0. Asymmetry only arises as a non-linear effect for short 
wavelength or large mass perturbations. 

These long-lived oscillations may help explain the presence of shells in elliptical galaxies. 
Quinn (1984) explained these patterns as transients that arise as an elliptical galaxy accretes a 
smaller system, and what we see is the track of a cannibalized dwarf galaxy, which forms shells 
through the phase wrapping process. Hernquist and Spergel (1992) suggested that these shells can 
also form through the merger of two equal mass spiral systems. The observed concentric structure 
of shells implies a large amount of dynamical friction and fast phase mixing, thus these shells are 
relatively short lived. Schweizer and Seitzer (1988) find that a large fraction of all ellipticals have 
shells, which implies a high accretion or merger rate with corresponding cosmological implications. 
From (|3|) we see that a sufficiently long wavelength perturbation does not exhibit dynamical 
friction and is subject to perturbative treatment. The diffuse remnant may continue to orbit 
for several dynamical times with little diffusion or damping. Furthermore, the periodic motion 
of particles in a shell can excite similar modes in the predator galaxy, which are also long lived. 
Therefore, a numerical simulation must necessarily take these excited modes into account in order 
to model such phenomena accurately. 

Under certain assumptions. Case and others have proven that the spectrum of eigenvalues 
is continuous except for a finite number of discrete points. For an Antonov stable system we 
know that they cannot contain positive imaginary components. The matrix analysis suggests 
that complex eigenvectors come in complex conjugate pairs, so we expect all modes to be stable. 
Landau damping then only occurs through the loss of coherence in a superposition of states when 
projected onto the density axis. But any periodic process would preferentially excite coherent 
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modes, as would be the case in satellites orbiting about larger galaxies. These modes can survive 
even after the satellite is disrupted or accreted. 

Landau's approach using Laplace transforms is an equivalent method to solve the initial value 
problem. Stability here implies the lack of frequencies with positive imaginary components. In 
the Jeans' argument, the damping occurs despite the existence of coherent van Kampen modes, 
since these are singular and thus violate the assumptions of the Landau analysis. In this paper we 
suggested that van Kampen modes need not in general be singular. This implies that there may 
exist poles in the inverse transform on the real axis, allowing undamped modes to exist. 

The success of the series expansion for the Jeans instability suggests that this scheme is a 
practical algorithm to calculate normal modes in arbitrary potentials. A systematic search for 
observable modes is now possible. A large density of states near certain eigenvalues might lead to 
easily excitable modes. 

A numerical eigenmode expansion of the full six dimensional system should be feasible if we 
choose the appropriate basis. Since our examples all have a basically Gaussian velocity space 
structure, together with spherical symmetry and power law radial dependence, the natural basis 
should involve Hermite functions in u-space, spherical harmonics for angular dependence and 
Bessel functions for the radial piece. One can envision a p-mode analysis of galaxies similar to 
their very successful application to stars. Observations of both the surface density and spatially 
resolved velocity perturbations can in principle supply us with a three-parameter data set, which 
we can compare to numerical calculations. This should enable us to determine the low order 
harmonics. Given infinite resolution, a three parameter observation allows us to infer the full three 
dimensional structure of the gravitational potential. From @ one expects that long wavelength 
perturbations should be longest lived, since they are the most linear and less subject to dynamical 
friction. This simplifies life for both the observer and the theorist, since a moderate resolution will 
pin down the fundamental modes and low harmonics. 



6. Conclusions 

The picture presented in this article gives a simple physical interpretation of perturbations 
about collisionless systems in terms of eigenmodes. I have presented some examples of periodic 
and stable oscillations in one dimensional collisionless systems, and discussed general features of 
these modes, which one can try to find in elliptical galaxies. 

We have explored a little of the structure of perturbations in velocity space, and saw 
that it lends itself to simple analysis when expanded in a suitable way. The analytic solutions 
are limited, but numerical studies can provide a detailed description of fundamental modes. 
An expansion in Hermite polynomials allows an explicit transformation of the perturbation 
equations into self-adjoint form for general integrable functions, which relaxes the standard 
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differentiability requirement. This argues that general stationary systems should also exhibit 
periodic nondispersive modes. 

The linear perturbation solutions of the six dimensional phase space are vast, and we have 
seen but a tiny sampling of its rich solution space. 
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